Numerical Integration Scheme Using Singular Perturbation Method

Some multi degree-of-freedom dynamical systems exhibit a response that contain fast and slow variables. An example of such systems is a multibody system with rigid and deformable bodies. Standard numerical integration of the resultant equations of motion must adjust the time step according to the frequency of the fastest variable. As a result, the computation time is sacrificed. The singular perturbation method is an analysis technique to deal with the interaction of slow and fast variables. In this study, a numerical integration scheme using the singular perturbation method is discussed, its absolute stability condition is derived, and its order of accuracy is investigated. INTRODUCTION Quite often, the solution of a state equation has some variables evolving in time faster than other variables, leading to the classification of variables as slow and fast. Such systems are called highly oscillatory system and the computation of their solutions is referred to as a multiple-timescale problem [1]. Standard numerical methods require a short step size to capture the dynamics of the fast variables. The computational cost of solving the entire system is dictated by the time scale of the fast variables and, hence, the numerical efficiency can become an important issue. Multirate methods were proposed to improve the numerical efficiency when solving highly oscillatory dynamical systems. The methods exploit the different time scales by using different step sizes for the subsystems. This kind of approach was applied for electric circuit simulation [2], molecular dynamics simulation [3, 4], and stellar problems [5]. There have been many attempts to apply multirate method to mechanical systems, especially multibody systems. In [6], partitioned Runge-Kutta method was employed to simulate the dynamics of vehicle systems that contain subsystems with high frequency response characteristics. The multirate method based on the Backward Differentiation Formula (BDF) was proposed and applied for the simulation of the aero-elastic model of helicopter [7]. In [8], the simulation of pantograph and catenary was conducted with a multirate method. In this study, a numerical integration method that utilizes the local linearization method [9, 10] and the singular perturbation method [11] to deal with the coupling between the fast and slow variables is introduced. The singular perturbation method is an analysis technique to deal with the interaction of slow and fast variables. The proposed numerical integration method can capture the fast dynamics using the local linearization method while the dynamics of the slow variable is computed by a conventional numerical method with the help of the invariant manifold of singular perturbation theory. The local linearization method is an exponential method which is based on the piecewise linear approximation of the state equation through a first-order Taylor expansion at each time step. The solution at the next time step is determined by the analytic solution of the approximated linear system. This paper discusses the absolute 1 Copyright c © 2013 by ASME stability condition and accuracy of the proposed numerical integration method and demonstrates its advantage by numerical experiments. SINGULAR PERTURBATION The singular perturbation model associated with a dynamical system is a state model where the derivatives of some of the states are multiplied by a small positive parameter ε [1]: ẋ = f (t,x,z,ε) (1a) ε ż = g(t,x,z,ε) (1b) It is assumed that the functions f and g are continuously differentiable and x ∈ Rn,z ∈ Rm. By setting ε = 0, Eqn. (1b) becomes 0 = g(t,x,z,0) (2) If Eqn. (2) has k ≥ 1 isolated real roots z = hi(t,x), i = 1,2, . . . ,k (3) then Eqns. (1a)-(1b) are in standard form and they reduce to ẋ = f (t,x,h(t,x),0) (4) This model is called the slow, reduced, or quasi-steady-state model. With the new time variable τ = (t− t0)/ε , the so-called boundary-layer model is defined as dy dτ = g(t,x,y+h(t,x),0), where y = z−h(t,x) (5) A geometric view of the singular perturbed system and the reduced model can be obtained by the concept of invariant manifolds. For simplicity, we consider the autonomous singularly perturbed system ẋ = f (x,z) (6a) ε ż = g(x,z) (6b) Let z = h(x) be an isolated root of 0 = g(x,z). Then, the equation z = h(x) is an invariant manifold for the system ẋ = f (x,z) (7a) 0 = g(x,z) (7b) When ε = 0, any trajectory starting in the manifold z = h(x) will remain in that manifold for all positive time. The dynamics in this manifold can be described by the reduced model as